%% create discretized maps and graphs of road networ
tic
clear
close all
clc

% set(0,'DefaultFigureWindowStyle','docked','DefaultFigureVisible','on')

% -----------------------
% DEFINE FOLDER LOCATIONS
folders;

%% code control

    % which cities to keep?
    popthreshold = 0.1;       % only load places with population above threshold
    pick_largest = 0.6;       % =1 to use most populated place in a cell to locate the node in the cell, otherwise choose centroid
    max_N_cities = inf;       % keep the max_N_cities largest cities
    N_cat = 10;                % for plotting only. number of population quantiles is N_pop_cat+1.

    % discretization
    % these parameters control the discretization of the road network
    diagonals = 1;    % include diagonals 
    x_ver_hor = 0;  % tolerance for vertical-horizonal links (higher=fewer links), from 0 to 1
    x_diag = 0;
            
    % transform roads into graph?
    road_network_graph = 1;  % to transform roads into graphs
    small_grid = 1; % indicator for Peru, Chile to look at the smaller grid 

    % add unconnected nodes?
    connect_all_nodes = 1;   % number of times the original network is made denser (times unconnected nodes are connected)
    largest_connected = 1;
    tol = 0.001; 
    
    % testing income data 
    winsorize = 1;
    
    % when to build for forests
    build_thresh = 80;
    
    % Block off the Amazon?
    brazil_forest = 0;
    
    default_cellsize = 0.5;   % default size is 50km by 50km (0.5)
    adjust_cellsize = 0;
    
    % which pop source to use?
    pop_source = 0; % 0 - default, SEDAC. 1 - WorldPop, 2 - Landscan
    
    % stock all code control variables
    code_control.build_thresh = build_thresh;
    code_control.datafolder_ne = datafolder_ne;
    code_control.popthreshold = popthreshold;       
    code_control.pick_largest = pick_largest;        
    code_control.max_N_cities = max_N_cities;    
    code_control.N_cat = N_cat;               
    code_control.road_network_graph = road_network_graph;
    code_control.connect_all_nodes = connect_all_nodes;   
    code_control.diagonals = diagonals;
    code_control.winsorize = winsorize;
    code_control.brazil_forest = brazil_forest;

%% load country list
country_list_short;


%% loop through selected countries
for country_n= 1:Ncountries
 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% load international boundaries %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    country_icc = char( countries(country_n) );  % country ICC code
    country = icc2name( country_icc );           % country name
    
    % get threshold from 
    x_ver_hor = get_threshold( country );
    x_diag  = x_ver_hor;
    
    % assert(x_diag>0)
    
    % add parameters to code control 
    code_control.x_ver_hor = x_ver_hor;  
    code_control.x_diag = x_diag;
    
    if ~ispc()   %% we are in a mac

        % folder to save the outcome
        datafolder_LAC_country =  [ datafolder_LAC,country_icc,'/'];

    else   %% we are in a pc   

        % folder to save the outcome
        datafolder_LAC_country =  [ datafolder_LAC,country_icc,'\' ];

    end
    
    % load country boundaries
    disp( ['Country: ',country] )
    file  = [ datafolder_LAC_country,'country_bounds.mat'];
    load(file)
           
%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% load populated places %%%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % load Population Data
   file  = [ datafolder_LAC_country,'popplaces_sedac.mat'];
   temp = load(file);   
   places_sedac = temp.places_sedac;
   clear temp
     
   % total population
   totpop = sum( cell2mat( {places_sedac.population} ) );
   
    % World Pop
    file  = [ datafolder_LAC_country,'popplaces_worldpop.mat'];
    load( file,'places_WorldPop' );

    % LandScan
    file  = [ datafolder_LAC_country,'popplaces_landscan.mat'];
    load( file,'places_LandScan' );


    
%%  %%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% bring income data %%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    
    file  = [ datafolder_LAC_country,'places_gecon.mat'];
    load( file,'places_gecon' );
    
%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% bring altitude data %%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    file  = [ datafolder_LAC_country,'relief_etopo.mat'];
    load( file,'relief_etopo' );     
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% load in travel speeds %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    file  = [ datafolder_LAC_country,'relief_etopo.mat'];
   
%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% discretize the map %%%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [ code_control,places_grid,gridmap,edges,places_gecon ] = ...
    brazil_grid( code_control,country_bounds,places_sedac,relief_etopo,places_gecon);

    % build everywhere - we don't impose that we can't build in 
    for j=1:length(gridmap)
        gridmap(j).build = 1;
    end
    
    % load in speeds data
    file  = [ datafolder_OSM, country_icc, 'coordinates_wtimes.csv'];
    osm_speeds = readtable(file);
    speeds  = cell2mat({osm_speeds.speed})';
    distances  = cell2mat({osm_speeds.travel_distance})';
    
    % Missing links that are neighbors can be walked (at very slow speeds)
    speeds(speeds<4) = 4;
    speeds = speeds';
    
    % get unique nodes from the gridmap 
    unique_nodes = ones( length(gridmap), 2);
    for n =1:length(unique_nodes)
            unique_nodes(n, 1) =gridmap(n).coordinates(1);
            unique_nodes(n, 2) =gridmap(n).coordinates(2);
    end
   
    unique_nodes_nat = unique_nodes;

    % move cities and discretize roads     
    places2 = move_cities( gridmap,places_grid,unique_nodes,unique_nodes_nat);  
    
    % no diagonals here since it is not a true grid 
    if strcmp(country, 'Brazil')
       code_control.diagonals = 0; 
    end
    
    [ discretized_roads,graph_export,unique_edges, dist_info ] = discretize_roads_speed_brazil( code_control,gridmap,places2, speeds, distances,edges,unique_nodes, country );        
    
    country_graph.places_grid=places_grid;
    country_graph.country_bounds=country_bounds;
    country_graph.gridmap=gridmap;
    country_graph.places2=places2;
    country_graph.discretized_roads = discretized_roads;
    country_graph.graph_export = graph_export;
    country_graph.edges = edges;  
    country_graph.unique_edges = unique_edges;      
    country_graph.code_control = code_control;
    country_graph.unique_nodes = unique_nodes; 
    country_graph.dist_info = dist_info;
    
  save( [ path_save_grids,country,'_grid_',...
    num2str( x_diag ),'_',...
    num2str( x_ver_hor ), '_',...
    num2str( default_cellsize ),'_connected', ...
    num2str( largest_connected ), ...
    '_LAC.mat' ],'country_graph' );   
    
end



        